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Fully doval^pod turbulant channal flow haa boon alnulatad nunarlcally at 
Raynolda nunbar 13800, based on centerline velocity and channel half width* 
The large-scale flow field has been obtained by directly Integrating the fil- 
tered, three-dlnenalonal, tine-dependent, Navler-Stokes equations* The snail- 
scale field notions were slaulated through an eddy viscosity nodal* The 
calcttlatlone were carried out on the ILLIiC IV conputer with up to 516,096 
grid points* 

The conputed flow field was used to study the statistical properties of 
the flow as well as its tine-dependent features* The agreenent of the 
conputed nean velocity profile , turbulence statistics , and detailed flow 
structures with experinental data is good. The resolvable portion of the 
statistical correlations appearing in the Reynolds stress equations are 
calculated* Particular attention is given to the exanlnatlon of the flow 
structure in the vicinity of the wall* 

1 * I ntroduction 

Large-eddy sinulation (LBS) is a relatively new approach to the calcula- 
tion of turbulent flows* The basic idea stens fron two experinental observa- 
tions* First, the large-scale structure of turbulent flows varies greatly 
fron flow to flow (e*g*. Jets vs. boundary layers) and consequently is diffi- 
cult, if not inpossible , to nodel in s general way* Second, the snail-scale 
turbulence structures are nearly Isotropic, very universal in character 
(Chapnan, 1979) and hence nuch nore anenable to general nodeling* In 1£8, 
one actually calculates the large-scale notions in a tine-dependent, three- 

*Portions of this irork were carried out while the authors held NRC 
Research Associateships at Ames Research Center. 


dinenalonal computatioa, using for ths Isrgs^scals field dynsmlcsl equsclons 
that Incorporate sluiple aodela for small'- scale turbulence. Only Che part of 
the turbulence field with scales that are small relative to overall dimensions 
of the flow field Is modeled. This Is in contrast to phenomenological turbu- 
lence modeling, In which all the deviations from the mean velocity profile are 
modeled. 

A typical LES calculation for wall-bounded turbulent flows Imposes a 
great demand on computer speed and memory. At present, therefore, the use of 
LES for practical engineering applications Is admittedly uneconomical. How- 
ever, for simple flows, such calculations are just within reach of the largest 
present computers. The information generated by these computations can In 
turn be used as a powerful research tool In studies of the structure and 
dynamics of turbulence. In addition, the various correlations that can be 
obtained from the computed large-scale flow field may be used ii.' developing 
phenomenological turbulence models for complex flows. These are the consid- 
erations that motivate the present development of the LCS method. 

The first application of LES was made by Deardorff (1970), whu simulated 
a turbulent channel flow at an indefinitely large Reynolds number. In this 
pioneering work he showed that three-dimensional computation of turbulence (at 
least for simple flows) is feasible. Using only 6,720 grid points, he was 
able to predict several features of turbulent channel flow with a fair amount 
of success. Of particular significance was the demonstration of the potential 
of LES for use in basic studies of turbulence. 

Following Deardorff 's work, Schumann (1973, 1973) also calculated turbu- 
lent channel flow and extended the method to cylindrical geometries (annuli) . 
He used up to ten times more grid points (65,336) than Deardorff and an 
improved subgrid scale (SGS) model. In addition to dividing SGS stresses into 
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a locally Isotropic part and an inhonc<teneoua parr , tie employed a separate 
partial differential equation for SGS turbulent kinetic energy. However, the 
added differencial equation did not Improve the results over the calculations 
in which only an eddy viscosity model was used (Schumann, 1975). 

GrStzbach and Schumann (1977) extended their channel flow calculations to 
account for temperature fluctuations and heat transfer. Later extensions by 
GrStzbach include calculations of secondary flows in partly roughened chan- 
nels, Inclusion of buoyancy effects, and liquid metal flows in plane channels 
and annuli. A recent review of this group's work in LES was given by Schumann 
et al. (1979). 

In all of the above computations, the dynamics of the inner region of the 
boundary layer (viscous sublayer and buiter layer) was essentially ignored. 
It is in this region that virtually all of the production of turbulence 
kinetic energy takes place (Townsend, 1956; Kim et al., 1971). Artificial 
boundary conditions in the logarithmic region wre used to simulate the Inner 
layers. Aside from the fact that these bojndary conditions are designed to be 
consistent in the mean with the law of tht wall, there Is little justification 
or experimental evidence to warrant their use for Che detailed flow field. 
However, the computations of Deardorff (1970) and especially those of the 
Karlsruhe group have produced successful comparisons with experimental data In 
the regions away from the walls. With a relatively modest number of grid 
points, they have been able to extract considerable information of practical 
Interest from their computations. 

The first numerical simulation of turbulent channel flow that computed 
rather than modeled the flow in the immediate neighborhood of the wall was 
that of Moln et al. (1976). In this calculation only 16 uniformly spaced grid 
points were used in each of the streamwlse, x, and spanwlse, z, directions 
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and 63 non-unlformly spaced grid points in the direction normal to the walls. 
The computational grid resolution in the lateral directions was Inadequate 
for resolving the experimentally observed coherent structures in the viscous 
sublayer. Nevertheless » the computations did display some of the well- 
established features of the flow in the wall region. The results of this 
computation were encouraging enough to Justify the undertaking of the present 
calculations . 

In this paper, we describe our numerical studies of Incompressible turbu- 
lent channel flow with up to 316,096 grid points. Particular attention is 
given to the investigation of the detailed flow structures. The Reynolds 
number, Re^ , based on shear velocity, u^ , and channel half-width was set 
at 640. The corresponding Reynolds number based on centerline velocity and 

channel half-width is about 13600 (Hussain and Reynolds, 1975). The results 
of this work can be summarized briefly by stating that, in the present compu- 
tations, the calculated mean velocity profile and turbulence statistics are in 
^ood agreement with the experimental data. The detailed time-dependent flow 
structures are strikingly similar to those observed experimentally. In 
addition, the resolvable portion of several statistical correlations which 
play an important role in phenomenological turbulence modeling are computed. 
These results tend to Indicate that the LES method can be used very effec- 
tively in supplementing laboratory measurements of turbulent shear flows. 

2 . Governing Equations for the Large-Scale Field 

In LES, each flow variable f is decomposed as follows: 

f - f -f f (2.1) 

where f is the large-scale component and f is the residual field. 
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Following Leonard (1974), we define the large-ecale field as 
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( 2 . 2 ) 


Where Is the filter function In the 1-direction and the integral Is ex- 
tended over the whole flow field. In planes perallel to the walls in which 
the flow Is statistically homogeneous, we use the Gaussian filter function: 

I /2 

Gj^(x^,xp - exp j^-6(x^ - , 1-1,3 (2.3) 

here, - 2h^ (Kwak et al., 1975), h^ is the computational mesh size In 
the 1-dlrectlon, and subscripts 1 and 3 refer to the streamwlse and spanwise 
directions, respectively. The corresponding integrals in (2.2) are extended 
over the entire (x^^iX^) plane. The width of the Gaussian function char- 
acterizes the size of the smallest eddies in the homogeneous directions that 
are retained In the filtered field (the largest eddies in the residual 


field) . 

Due to variation of turbulence length scale in the direction normal to 
the walls, X2 , one should use a filter with a variable width, A2(x2) . In 
this direction a sectlonally continuous “top hat" filter function was used. 
Let X2^ be the location of the computational grid point In the ver- 

tical direction; we define tJie filter function 63 for the control volume 


surrounding X2^ as follows 


(A'^'Cx^) A~(x 2))"^ for X2 “ a“’(x 2 ) < x^ < *2 + a'*‘(x 2) 


G2(x2»xp - 


U.'*) 


for x^ > x^ A (x^) and x^ < x^ ~ A (X2) 


where 
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The functioas and tT are aectlonally conscant functions of X 2 > 
therefore, in the open neighborhood surrounding each computational grid point, 
X 2 - A“ < X 2 < *2 + dA‘^/dX 2 ■ dh~/dX 2 "0. An important consequence 
of this property of and A~ and the form of and G-^ (function of 
Xi - x^, i ■ 1,3) is the commutivity of the filtering operation and partial 
differentiation operators in these neighborhoods and in particular at the 
computational grid points (see Hoin et al. (1978)), i.e., 


IT 31 

11^ “ Iiq- 


(2.5) 


Note that, with the appli<^ation of G 2 » the filtered variable f will be 
sectionally continuous and the filtering in (2.3) is interpreted as an average 
over grid spacings in the X 2 -directlon. This is a step prior to complete 
discretization of flow variables for numerical computation (Section 3) . 

Schumann (1973) and co-workers use a filtering function similar to (2.4) 
in ail spatial directions. When applying this averaging process to the 
Navier-Stokes equations, they evaluate the integral in the direction of the 
derivatives analytically and then have to deal with averages over the faces of 
the control volumes. This process introduces three types of surface-averaged 
as well as volume-averaged variables. The extra variables have to be related 
to each other in some way. 

Now, applying the filtering operation (2.2) to the Incompressible Navier- 
Stokes and continuity equations, we obtain the dynamical equations of the 
large-scale flow field. 
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( 2 . 6 ) 
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where we have decompoaed as in (2.1) and 
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here, Che variables are non-dioenslonalized using Che channel half-wldch 6 
and Che wail shear velocicy u^. The calculaclons were carried ouC for a 
fixed screamwlse mean-pressure gradienC which is accounCed for by Che 
Cenn in Che mor:enCun equacion (2.6). 

There are cwo polncs associaced wich Eq . (2.6) ChaC require furCher 
explanaClon. FirsC, Che convecCive cerm, *3^ wriccen in che 

J 

equlvalenC buC more cumbersome form 


■ ‘tJk“A * 3^ 7 “j“j 

This was done because ic can be shoim (Hansour eC al., 1977) ChaCi wich chis 
form in Che absence of clme-dif ferencing errors, conservacion of energy, 
momenCum, and circulacion in inviscid flows will be obcained when vircuaily 
any difference scheme is ftpplled co (2.6). Second, Ic should be noced chaC 
Che so-called Leonard (1974) sCress cerm 
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(2.9) 



U U u u 

1 J 1 J 


is not equal to zero. One has the option of calculating the tens with double 
bars explicitly, or, as Deardorff (1970) has done, to incorporate in 
whatever aodeling assumption is used for (see the next section). In the 
present work, with respect to the Gaussian filter in the horizontal directions 
where the partial derivativea are calculated pseudospectrally (bee. S), we 
have chosen the former option. In Sec. 6 we shall show that X^j can be 
quite significant^ hence. Including It with Is not recommended. With 
respect to the top**hat filter In the X 2 -direction where the derivatives are 
evaluated by second-order finite-difference schemes, the latter option was 
chosen. Here, It can be shown that (Shaanan et al., 1973) X^j is of the 
same form and order as the truncation error of the finite-difference scheme; 
hence Its explicit calculation la not Justified. However, when higher-order 
finite-difference schemes or spectral methods are used to evaluate the deriv- 
atives In the X 2 ~directlon, X^j should be calc\ilated explicitly. 


3 . Residual Stress Model 

The basic Idea behind large-eddy simulation Is that the large-scale 
motions, which are calculated explicitly, provide most of the important 
turbulent transport and, hence, the Influence of the small eddies can be 
modeled relatively crudely. In the present calculations, we have used an eddy 
viscosity model for similar to that used by Schumann (1975). 



(3.1) 
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and < > denotes the average over a plane parallel to the walla. The small* 

scale eddy viscucttles and represent the action of the unresriv&l 

scales of motion on those chat are resolved. Hence, as the resolution gets 
better, and should get smaller. Since the filter widths represent 

the largest length scales that are not resolved, these widths are the charac* 
teriatic size of the largest (and hence moat important) residual motions. 

The first term in (3.1) with 

- (C A)^ /2(S >HS >) (3.2) 

T s IJ IJ IJ IJ 

is Smagorinsky's model and can be derived from equating tne subgrid scale 
(SGS) production and dissipation in homogeneous turbulence. This model was 
used successfully in the numerical simulation of the decay of isotropic 
turbulence by Mansour et al. (1977) and by Oeardorff (1970) (with S^j - 
replaced by in the calculation of tlie cere region of turbulent channel 

flow. In Che expression for Vj (Eq. (3*2)), A is the characteristic 

length scale of the largest subgrid scale eddies, here assumed to be (Uear- 
dorff, 1970) 

A - (A^A2Aj)^/^ (3.3) 

is a dimensionless constant, and A^ 4s the filter width in the i- 
direction. In addition, in order to account for low Seynolds number SGS 
turbulence near the wall, the above expression for A was multiplied by the 
Van Driest (1936) exponential damping function, 1 - exp(-y'^/A^) , with - 
25 and • y u^/v, the distance to the nearest wall in the wall units. In 
all the calculations reported here which were performed with different grid 
sizes, the value of ■ 0.065 was used. Numerical experiments indicated 
that the use of a value much larger than this one resulted in excessive 
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damping of Che resolvable turboleoce* When lower voluea of Cg %iere uaed, 
exceaa energy acciLiulated near the high wave number end of one^dimenalonal 
energy spectra. In general » however, the compuced turbulence Intenalcles were 
rather Inaenaltlve to small perturbations 201) of C^. Note that the 
above value of Cg corresponds to Cg ■ 0.1, uaed by Oeardorff (1970), if 
in (3.3) is replaced by h^« 

Nuar Che wall the important large scale atructures are the "streaks" 
(Kline et al., 1967). These structures are relatively finely spaced in the 
spanwise direction. Their mean spacing characterizes the length scale of the 
eddies in the viscous sublayer (and hence the thickness of the viscous sub- 
layer) . Thus in a calculation with inadequate resolution in the spanwise 
direction, the thickness of the viscous sublayer will probably be larger than 
its physical counterpart. This in turn will lead to lower gradients of the 
computed mean velocity profile and consequently insufficient production of the 
resolvable turbulent kinetic energy. Therefore, in order to account for the 
effect of some of the streaks that reside in SGS motion on the mean velocity 
profile, the second term In (3.1) was introduced in the model for t^j. \s 
was mentioned earlier, Schumann (1975) ha.) also decomposed SGS stresses Into 
two parts, as in (3.1). The first was to account for locally Isotropic SGS 
stresses, and the second to account for inhomog^neities due to the nonzero 
component of mean strain. 

In the present study, the eddy viscosity, Vj, is defined as follows; 

V* - cmy n < S,' > < S, > (3.4) 

T 3 IJ IJ 

where c ■ 0.065 is a dlmenalonl<:'i«.' constant and D ■ 1 - exp(-y'*’^/A^‘^) is a 
VMll-damping function with « 25, as before. It should be pointed out 

that the characteristic length scale associated with is ^j, the filter 
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width In the spanwlse direction. At the reeolution in the s-direction is 
improved, v!^ will epproech eero. Moreover, it ehould be noted that, due to 
its functional form (function of y only), Vj doee not appear directly in 
the governing equationa for the reaolvable portions of turbulence stresae.s and 
hence doee not contribute to the dissipation terms in these equationa. This 
is in contrast to v^, which will suppleuent the molecular viscosity as a 

dissipating agent for the resolvable turbulence stresses. However, con* 

1 2 

tributes to the dissipation of mean kinetic energy, u > , and therefore 
indirectly to the production of resolvable turbulence stresses. The value of 
c in £q . (3 .A) was chosen to be 0.063 from a numerical experiment. It is 
approximately the minimum value with which the resolvable turbulent kinetic 
energy can be maintained (i.e., it did not decay indefinitely). This numeri- 
cal experiment was performed with one of the computations reported in Table i 
(Case 1), but the same constant was used in all the other calculations repor- 
ted here. 


4 . The Computational Grid 

Three faccors Influence the choice of the computational grid. First, the 
mesh size should be small enough to resolve tfte important scales of motion la 
the flow. Secor.j, the computational dooialn should be large enough tliat arti- 
ficialities of the boundary conditions do not influence the statistics of the 
solution in an undesirable way. Third, the availability of computer resource 
restricts the size of calculation that can be done. 

In the direction normal to the walls (-1 ^ y ^ 1)^ 63 grid points with 
non-uniform spacings were distributed. The following transformation gives the 
location of grid points in this direction: 

^For notatlonai simplicity, we occasionally shift from (^l,X 2 ,x^), 
(Uj^,U 2 ,U 3 > to (x,y,z), (u,v,w). 



( 4 . 1 ) 


Xj ■ y taah tanh ^(a) 

where 

Cj - - 1 + 2(J - D/CNj - 1) . j - ^2 

N2 Is the total number of grid polnta In the y-dlrectlon* Here a la the 
adjustable parameter of the transformation (0 < a < 1); a large value of 
a distributes more points near the walls. In our computations we have used 
a “ 0.98346, N2 * 63. This value of a was selected so that the above grid 
distribution In the y-dlrectlon Is sufficient to resolve the viscous sub- 
layer (y"*" < 5) . 

The selection of the length of the computational box In the streamwlse 
and spanwlse directions Is Initially guided by the two-point correlation 
measurements of Comte-Bellot (1963). Her data show that the correlation 
between velocity fluctuations at two points away from the walls* separated 
In becomes negligible at an x^ separation of 3.26. The correlation 

between motions at two points (away from the walls) separated in X3 ■ becomes 
negligible beyond an x^ separation of 1.66. thus, if we wish to employ 
periodic boundary conditions in x^ and x^ directions, we must choose a 
computational domain approximately twice as large as these dimensions. This 
is to prevent these simple but artificial boundary conditions from seriously 
influencing the results (Schumann, 1973). It should be noted, however, that 
the computed two-point correlation functions provide sufficient information 
regarding the adequacy of Che computational domain. If, for example, in the 
x^^-direction the length of the computational box, , is too short, the 
computed profile(s) of K^^(y,r) does not decay sufficiently in the neigh- 
borhood of r • Lj^/2, and hence should be Increased. 

Near the wall data are not available . 


] 
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In the vail region, the Important large-acale atructures are the 
“atreaka” (Kline et al>, 1967). Theae have a mean apanvlee apaclng corre- 
apondlng to » 100 , vith the moat probable apanvlae aparxng, A^p, about 
80 in the vail unlta. In addition, Kline et al. (1967) and Clark and Markland 
(1970) occaalonally obaerved U-ahaped vortlcea in the inner region. In the 
atudies of Clark and Markland, the average atreamvlae apaclng of these 
structures vas found to be A^n^ > 440. For the present computation at Ke^ ^ 
640, these correspond to dimenslonlesa apacinga of 

^Im “ 440/640 - 0.687 

^ 3 p ■ ®0/640 - 0.125 

Table 1 shows the characteristics of the computational grid networks in 
four different calculations reported here. In this table, is the number 

of grid points; is the length of the computational box, h^ is the grid 

spacing in the i-direction, and subscripts 1 and 3 refer to the streamwise 
and spanwise directions, respectively. The ^/6 and i* 2 /»S entries in Table 
1 show tijat, except for case 1, where Lj < 3.26^, the size of the computa- 
tional domain in all other cases appeals to be large enough to accommodate the 
important large eadies. Furthermore, since the pseudospectrai method (Section 
5 ) is used to approximate the derivatives in the streamwise and spanwise di- 
rections, the computational grid resolution (at least for Cases 2, 3, and 4) 
appears to be Just adequate to resolve structures with A|^|^ and A^p spacing 
in the Xj^- and x^-directions , respectively. It is emphasized that the ahuve 
values for A|^|^ and A^p are based on an ensemble of measurements, and, at a 
given instant, structures with a finer spacing than A^p and Aj^^^ can be 

^In this case, the computed two-point correlation functions R 33 (y>rj) 
(for y > .26) Indicate that L 3 is not sufficiently long. 
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formed which cannot be resolved with Che current grid resolution. Thus, we 
cannot expect the present calculations to reveal the streaky structures In the 
viscous sublayer with mean spacing equal to As we shall see, however, 
calculations do produce screaks at the finest scale permitted by the grid. 

Finally, we mention that the grid meshes for pressure do not coincide 
with grid meshes for velocities. Grid points for P are located midway between 
those for u^^. The equation of continuity Is enforced at node points for 
P, whereas the momentun equations are evaluated at node points for u^^. Note 
that, In contrast to the conventional staggered grid system (Harlow and Welch, 
1963), In which the three velocity components are defined at different node 
points, In the present grid system all the velocities are defined at the same 
grid points. This will allow for convenient application of the wall boundary 
conditions • 

5. The Numerical Method 

Partial derivatives In the and x^ directions were evaluated pseu- 
dospectrally (Orszag, 1972). This Involves Caking the xj^ (or X 3 ) t'ocrier 
transform of the function to be differentiated, multiplying the result »y 
Ik^ (or Ik^) , where k^^ (k^) Is the wavenumber In the (X3) direction, 
followed by inverse transformation to get the desired derivative. This method 
lus the advantage that it handles the high wavenumber components of the func- 
tion precisely. Thus, the use of the pseudospectral method In the x^^ and 
Xj directions gives us the best possible resolution (with a given number of 
grid points) In these directions. Partial derivatives in the X 2 direction 
were approximated by central difference formulae. These will be described 
below. 

The time advancement was made using a semi- Implicit method (Moin et al., 
1978). The momentum equations (2.6) were recast In the form 
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(5.1) 
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where contains all of the terms in (2.6) that are not in (5.1). For 

discretization in time, we used the Adama-^Bashforth method for and the 

Crank-Nlcolson method for the remaining terms in the right-hand side of Eq . 
(3.1). For convenience, we evaluated < Vj > and at the old time step 

n. 


In Eq. (3.1), includes the term 
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The computation of this term can be accomplished by first calculating the term 
under the large overbar, taking the Fourier transform with respect to 
and x^i multiplying by the Fourier transform of the Gaussian filter func- 
tion, and then inverse transforming. 

Next, we Fourier transform the resulting equations in the x^ and X 3 
directions, this converts the set of partial differential equations (5.1) to 
the following system of ordinary differential equations, for every pair of 
Fourier modes k^ and k^ with y * X 2 as the Independent variable ; 
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Here, 6^(1 • 1,2,3) are known functions of Re^, < >”, and and 

A 

represent the terms involving pressure and velocity field at time 
steps n and n 1 • 

The following central difference formulae. 
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were used to approximate 
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respectively, in Eqs. (5.2). Here J denotes the velocity mesh point yj , 
q the pressure mesh point y^ , and " yj ” yj-l* resulting set of 
equations, together with the equation of continuity evaluated at the pressure 
node points , 


■'n-fl “n+1 



(5.4) 


leads to a system of algebraic equations for the Fourier transform of the de- 
pendent variables at the new time step. This system is of block- tridiagonal 
form and can be solved very efficiently (see below) . No-slip boundary con- 
ditions on velocity were used at the walls (y * * 1 ) ond periodic boundary 
conditions were incorporated in the Xj^ and X 3 directions. Note that 
pressure wall conditions are not necessary; only velocity boundary conditions 
are sufficient to close the system of equations (see Moin and Kim, 1980). 

In the present calculations, the core memory of tlie ILLIAC IV is large 
enough to hold only a few planes of the dependent variables. Therefore, it is 
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Important to manage efficiently the transfer of data between the core and disk 
memory, where the entire data base resides. A detailed description of the 
data-management technique used here Is given In Kim and Moln (1979). Here, we 
briefly outline the essential steps. At each time step, the system of alge- 
braic equations just described Is solved by two separate passes through the 
data base. In PASS I, the right-hand side of these equations Is computed. 
This Is accomplished by transferring two (x-z) planes of the Independent 
variables from the disk memory to the core memory to be processed by a double- 
buffer scheme. In this manner, all the (x-z) planes are transferred to the 
core, two planes at a time. 

In PASS 2, the block-trldlagonal system must be solved for each kj^ and 
k^. In this pass, (y~k^) planes of the right-hand side vector that were 
computed in PASS 1 are transferred to the core memory. Due to the limitation 
of the core size of the ILLIAC IV, a special algorittin had to be developed to 
solve the block-tridiugonal system of equations. For each and this 
algorithm requires 6'^6 N 2 floating-point arithmetic operations, in contrast 
to 376 N 2 operations for the conventional block-trldlagonal solver (Merriam, 
197ti). The extra computations are necessary In order to avoid the extra 1/0 
passes that would otherwise have been necessary. 

With a full use of the parallel processing capabilities of the ILLIAC IV 
computer and the above data-management technique , the computer time per time 
step (CPU and I/O time) was 22 sec for 63 x 64 x o4 grid-point calcula- 
tions and 36 sec for the computations with 63 x 64 x 128 grid points. For 
tlie calculations shown in Table 1, the dimensionless time step, At, was set 
at O.OCI for Case I, 0.00075 for Case 2, and 0.0005 for Cases 3 and 4. 
Throughout the computations , the value of 


C(t) 


Max 




+ 





I 

) 


never exceeded 0 .35 . 
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The initial condition for Case 1 in Table 1 was obtained by assigning the 
final velocity field described in Kim and Moln (1979) to the corresponding 
grid points used here. For Case 2, the final velocity field from Case 1 was 
simply assigned to the alternate mesh points. The values at the intermediate 
points were obtained by Fourier interpolation. A similar procedure was used 
to generate the initial velocity field for Cases 3 and A. 

6. Mean Velocity Profile and IVirbulence Statistics 

Starting from the initial velocity field, for each case, the governing 
equations were integrated forward in time until the numerical solutions 
reached statistically steady states. These equilibrium states were identified 
by approximate periodicity of the horizontally averaged turbulence stresses in 
time. Next, in order to obtain better statistical samples, the equations were 
further integrated in time and a running time average of the horizontally 
averaged turbulence quantities was calculated. For each case, the calcula- 
tions were considered to be complete when the time-averaged turbulence quanti- 
ties became stationary. The total time of integration and the averaging time 
for all the computations reported here are shown in Table 1. 

Figure 1 shows the mean velocity profile < IT > (unless otherwise stated 
in this section, < > indicates horizontal as well as time averaging) for 
ail cases reported in Table 1. The calculated mean velocity profiles are in 
good agreement with each other as well as with the experimental data of 
Hussain and Reynolds (1975). It appears that, at least for the different com- 
putational grid networks considered here, the law of the wail and the corres- 
ponding logarithmic layer and von KirmSn constant can be predicted with 
virtually no dependence on the grid resolution. This is particularly signifi- 
cant in light of the fact that the characteristic length scale of the subgrid 
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scale model used here Is a funccloo of computational grid resolution. Figure 
2 shows the profiles of the correlatloa coefficient between the resolvable 
streamwlse and vertical components of turbulence fluctuations , 


. -”2 .1/2 . -2 a/2, 

< u V >/I< u > < V > I 


wnere u ■ u ~ < u >. The calculated profiles for all the cases reported In 
Table I are In good agreement with each other and with the experimental data 
of Sabot and Comte-Bellot (1976). The results presented In Fig. 2 together 
with those in Fig. 1 establish some confidence in the reliability of the sub- 
grid scale model used in this study. 

In the remainder of this paper we shall present, in some detail, the 
results obtained from Case 4 in Table 1. The computational grid resolution 
for this case is better than for all the other cases; and, as will be shown 
later in this section, the computational domain appears to be large enough to 
Include the Important large eddies. 


6.1. Turbulence Stresses 

Vertical profiles of the resolvable mean Reynolds shear stress, < u v>, 

and the total Reynolds shear stress, < u v > -f < t ^2 shown in Fig. 

3. These profiles indicate that the average Reynolds shear stress profile has 

attained the equilibrium shape which balances the downstream mean pressure 

gradient in the regions away from the walls. In the vicinity of the walls, 

the viscous stresses are significant, and they, together with the total 

Reynolds stress, balance the mean pressure gradient. The symmetry of the 
•• 

|< u V > + 1 ^ 2 ! profile about the channel centerline Indicates that the 
total averaging time and statistical sample are adequate. Moreover, it should 
be noted that the subgrid-scale contribution to the total Reynolds stress is 
significant only in the vicinity of the walls. 
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Figures 4 and 5 show Che profiles of the dimensionless resolvable Cur- 
bulence Intensities. For comparison, some of the available experimental data 
over a range of Reynolds numbers are also shown. Once again, the symmetry of 
the calculated turbulence Intensities about the centerline of the channel 
Indicates that the total averaging time was sufficient for an adequate statis- 
tical sample. The overall agreement of the computed turbulence Intensities 
with the experimental measurements Is good. In Fig. 5, the resolvable tur- 
bulence Intensities In the vicinity of the lower wall are plotted vs. y*** * 
y^u^/v. In spite of large differences among various measurements, the maximum 
of the computed < u is located at a distance farther away from the 
wall (y^ ^ 30) than those of the measured turbulence Intensltlea (y*^ ^ 13- 
20). In addition. In the lomedlate neighborhood of the wall, an appreciable 
fraction of the vertical component of turbulence Intensity appears to reside 
in the subgrld-scale motions. It should be noted that. In contrast to the 
turbulence shear stress, in order to deduce the subgrid scale contribution to 
turbulence Intensities one has to obtain an estimate for the kinetic energy of 
SOS stress, and use It In Eq . (2.8a). Due to the lilgh degree of ani- 
sotropy In the channel flow, especially In the vicinity of the walls, we have 
been unable to obtain a reasonably accurate estimate for 


6.2. Two-Point Correlation Functions 
Two-point correlation functions 


Rli(y.ri) 


«• •• 

< u^(x,y,x) u^(x + rj^,y,z) > 

^ 

< u^ (x,y,z) > 


and 


< u.(x,y,z) u,(x,y,z r,) 


“ < u^'"(x,y,z) > 
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for 1 ” 1(2,3 (no summation) are plotted in Figs* 6 and 7 at four vertical 
locations. These profiles show that In general, for small separation dl8~ 
tances, the correlation for the velocity In the direction of the displacement 
Is larger than the corresponding transverse correlations. In addition, the 
longitudinal correlation In the streamvlse direction extends over much longer 
distances than do all other correlations. This result was also obtained by 
Deardorff (1970). 

The slow decay of " 0.025, r^) for Increasing r^ indicates that 

near the wall the eddies are highly elongated In the streamwise direction. On 
the other hand, the profiles of R^^Cy.rj) show that the spanwlse extent of 
turbulence structures near the wall is much smaller than for those away from 
the wall. Thus, it appears that, in accordance with the experimental observa- 
tions, near the walls the computed flow field consists of elongated streaky 
structures. The structure of the flow field will be examined in some detail 
in the next section. 

For comparison, in Figs. 6 and 7, the profiles of 
at y^/6 » 0.11, 0.44, and 1.0 from Comte-Bellot's (1963) measurements are 
included. Note that the computed and measured correlations are obtained at 
slightly different vertical locations. The correlations were calculated at 
selected points in the y-dlrectlon, and the comparison is made at the loca- 
tions where the y-coordlnate of the computed and measured correlations were 
closest to each other. For small values of the non-dimensional ized separation 
distance, r^, the measured correlations, *^n(yv*^l)» smaller than the 

computed ones, whereas for larger values of r^ the reverse is true. At 
small values of r^^ , the discrepancy between the computed and measured cor- 
relations is due to the fact that the measurements were made at a much 
larger Reynolds number (Re ■ 133,000) than in the present simulation (see 
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Batchelor, 1953) • However, Che cause of the difference between the computed 
and measured profiles for larger values of r^ is not clear* Comparison of 
Che measured profiles with the computed ones obtained from calculations with 
larger computational box lengths in the streamwise direction shows no improve- 
ment. Thus, the streamwise extent of the computational box does not appear to 
be a factor here. However, possible inadequate resolution in the x-direction 
may suppress the formation of some small scale structures. In this case, these 
structures could conceivably combine to form eddies that have long streamwise 
extent throughout the channel cross section. In addition, it should be men- 
tioned that Comte-Bellot's axial two-point correlation data were obtained by 
traversing one probe downstream of another probe. With this procedure, the 
measured Hj^^(y,rj^) may be contaminated with errors, due to the effect of the 
wake of the upstream probe. However, the probe interference effect should be 
significant only for small separation distances. 

In Fig. 7, the profiles of also compared with the mea- 
surements of Comte-Bellot . Aside from the Reynolds number effect for small 
values of r ^ , the agreement between the computed and measured correlations 
is good. Finally, the other two spanwise correlations, 

R|j 3 (y,r^), were also compared with the corresponding ones measured by Comte- 
Bellot. The measured R22(yi>^3) ^33(y»>^3) systematically lower 
than the computed ones for the values of r^ for which these correlations 
have appreclaoie magnitude. 
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6.3. Skcvnes* and Flatneat Factor* of Raaolvabl* TUrbulanca 


The velocity akewnesa aad flatneaa factor* which are defined as 

—**3 —"A 

< u /* > _ < u ^ > 

S(u ) - - j - H - y and F(u") - — j (1 - 1,2,3; no aummation) 

< ' 

respectively, are plotted in Fig. 8. The flatneaa fact'.ra of all the velocity 
conponents reach their maxima at the wall. This indicates that in the vicin- 
it; of the wall the turbulence is highly intermittent. Throughout an appreci-* 
able portion of the channel cross section, F(u^) and S(u^) are approximately 
equal to three and zero, respectively* These values correspond to the flat- 
ness and skewness factors of a Gaussian distribution. However, Kreplin and 
Eckelmann (1979) measured the U3 probability distribution and have shown 
chat it is not Gaussian, even though the values of S(u^) and F(U3) correspond 
to that of a Gaussian distribution. 

Near the wall, S(u^) is positive, whereas away from the wall it is 

negative. This indicates that near the wall the large-amplitude u fluc- 
tuations are primarily due to arrival of high-speed fluid from regions away 
from the wall. On the other hand, away from the wall the large-amplitude 
u fluctuations are most probably associated with low-speed fluid leaving the 

wall region. These observations are in agreement with the experimental find- 

•* 

ings of Brodkey et al. (1974) and with contour plots of IT 7 from numerical 
simulation of turbulent channel flow (Kim and Moin, 1979). however, as will 
become clear below, the precise vertical location (in wall units) of the 
crossover point in the present calculation is in disagreement with experimen- 
tal data. This discrepancy is probably due to inadequate grid resolution In 
the computations. 

In Fig. 8, the profiles of skewness and flatness factors from measure- 
ments of Comte-Bellot (1963) and Kreplin and Eckelmann (1979) are reproduced. 
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The overall agreement between computational and experimental data la good. 
This is particularly encouraging cuuaidering the significant contribution of 
small scale turbulence to these quantities and the difficulties associated 
with their measurements. 

6. A . Resolvable Vorticlty Fluctuations 

Figure 9 shows the profiles of the non-dimensional rms vorticlty fluc- 

2 1/2 

tuaclons. The spanwlse component of vorticlty fluctuations , ^ > , 

attains its maxlmian at the wall and decreases monotonically towards the chan- 
nel centerline. The profile of the rms streamwise vorticlty fluctuations, 
2 1/2 

< > , also attains its maximian at the wall but, in addition, displays a 

local maxlmu^i at y^ > 30. At about the same location, the peak of the rms 

vertical component of the vorticlty Is located. The mechanics underlying this 

behavior of the profile of < > ' will be discussed in tlte next section. 

It Is Interesting to note that, In spite of large differences between 

2 I '2 

different components of rms ^ ^ ' near the wall, away from the wail 

(y^ > 70) they are virtually identical. This is in contrast to rms velocity 
— "2 1/2 

fluctuations, < u^ > . The difference between the two may be explained by 

noting that the relative contribution of small scales to vorticlty fluctua- 
tions Is significantly larger than their contribution to velocity fluctua- 
tions, and away from the walls the small scales tend to be Isotropic. 
Ex'^loitlng the "Isotropy" of vorticlty fluctuations may be very useful in 
statistical analysis of c-urhulent shear flows. 

The limiting wall value of vorticlty fluctuations In the present calcula- 
tions are 0.20 and 0.13 tor < >1/2 ^ ^2 >1/2^ respec- 

tively. The agreement of these computed values with experimental measurements 
(see Kreplln and EcLelmann, 1979, for data from several measurements) is sat- 
isfactory . 
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6.5. Sf ti>tlct Involving RB>olv«bI« Pftuf Pluctuatlont 
The rooc*me«n equere velue of the reeolveble Mall preaaure fluctuations, 
< is 2.05 for both walls. This value is in fair agreement with 

the values of 2.64 obtained by Wlllaarth and Wooldridge (1962), 2.31 by Uiil- 

marth (1965), 2.6 by Elliot (1972), and of 2.0 and 2.5 reported by Corcos 

(i962) for fully developed pipe flows. However, the coisputed values are con- 
siderably lower than the 3.59 obtained by Blake (1970) or Eomerling's (1973) 
measurements (see Uillmarth (1975)).^ The measurements of Blake and Emerling 
were made with a pinhole microphone and by optical techniques, respectively. 
Therefore, in these experiments the smalier-scale pressure variations are 
expected to be better resolved. Thus, In view of these two experiments, it 
appears that an appreciable portion of the pressure fluctuations may reside in 
subgrid scale motions. 

The profiles of the diagonal elements of the resolvable pressure strain 
correlation tensor 



are shown in Figs. 10 and ll. These terms govern the exchange of energy 
between the three components of resolvable turbulence kinetic energy (hinze, 
1975). The negative sign for (no summation) indicates loss, or transfer 

of energy from < u^^ to other components, whereas a positive sign 

denotes energy gain. These profiles show that, except in the close vicinity 
of tile wall, as expected, the streamwise component of turbulent velocity 
fluctuations transfers energy to the cross stream components, however, very 
near the wall, there is a large transfer of energy from the vertical component 

*See Wlllmarth (1975) for the above references, as well cs values from 
other sources. 



of turbulence Intensity to the horizontal components. In this work, we shall 
refer to this phenomenon as the “splatting" effect. It will be shown that the 
splatting effect is an important property of the flow in the vicinity of the 
walls and should be taken into account in the modeling of near-wall turbu- 
lence. In fact, this phenomenon was noted in a previous study by Daly and 
Harlow (1970), who Included a term in their statistical model of <|)^j to 
account x its effect. 

The pt' xile of the off-diagonal element together with the pressure 

a 

diffusion term - < Pu > and their sum 

3y 


P 


uv 


- - < u > 


in the vicinity of the lower wall are shown in Fig. 12. The last term, 
appears in the governing equation for resolvable turbulence shear stress. 

g •• 

< u V >. The components of ♦ 12 ’ ~ Ty ^ ^ have comparable 

magnitudes and, as will be shown below, near the wall they provide important 

(t 

contributions to the governing equation for < vl'v >. As expected (Uinze, 
1975), except in the immediate neighborhood of the wall, the sign of 
opposite to that of < u v >. However, near the wall, where the splatting 
effect is present, «j>]^2 same sign as < u v >, thus contributing to 

the production of turbulence . 


6.6. Resolvable Turbulence Intensity and Shear Stress Balance 
In phenomenological turbulence modeling, the objective is to construct 
rational models for the correlations that appear in the governing equations 
for the Reynolds stresses. For the resolvable portion of the flow field in 
channel flow, these equations are: 
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It should be emphasized that the above equations are for the resolvable 
portion rather than the total turbulent stresses. However, considerable in- 
sight into the mechanics of energy transfer and the relative Importance of 
various terms in the Reynolds stress equations may be gained from the corres- 
ponding terms in these equations. 

In Figs. 13-17, all the correlations appearing in the above equations, 
and in the governing equation for the resolvable turbulent kinetic energy , 

, are plotted in the vicinity of the lower wall (y"^ 
< 90). The first term in the right-hand side of Eqs. (6.1) and (6.4) is the 
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production term. In Figs. 13-17, the remaining five terms in the right-hand 
side of each equation are labeled as convection, velocity-pressure gradient 
(VPG) , diffusion, cascade, and dissipation. The last term in each equation, 
A^j, la a relatively complicated expression involving v^. They would be 
identically equal to zero if were a constant. These terms were cal- 

culated and found to be negligibly small compared to the other terms in each 
equation. 

In Figs. 13 and 17, the production and dissipation are clearly the domi- 
nant terms in most of the region shown. In the lomedlate neighborhood of the 
wall, however, where the production term is small, viscous diffusion carries 
sufficient energy inward to balance the large viscous dissipation there. In 
addition, it can be seen that, aside from the close vicinity of the wall (y 
< 15), energy is convected from the wall region, where production is high, to 
Che regions away from the wall. 

The velocity-pressure gradient terms make large contributions to the 
balance of the governing equations for normal and spanwise components of tur- 
bulent kinetic energy. Near the wall, the triple correlation term (convec- 

—2 

tion) and pressure strain and pressure diffusion terms in the < U2 > 
equation are very significant. In particular, the reduction of the normal 
component of turbulent energy due to the splatting effect mentioned above is 
compensated by the pressure diffusion term. 

Near the wall (y^ < 25) , in the dynamical equation for Reynolds shear 
stress, < u V >, the triple correlation and velocity-pressure gradient, 
are the dominant terms. However, for y"*" > 25, the contribution of 

is small and the production term takes on a more active role. Moreover, 
it should be noted that in this equation the viscous diffusion and "dissipa- 
tion*' terms appear to be negligible. This in turn implies that the governing 

M 

equation for < u v > has a hyperbolic character. 


28 



In Eqs. (6*2)-(6«4) and in Figs* 13-17, the teraa Involving molecular and 
eddy viscosity were combined* In Fig. 18, the dissipation of turbulent kin- 
etic energy due to molecular and eddy viscosity are plotted separately. It 
can be seen that, near the wall, the dissipation due to eddy viscosity is neg- 
ligible compared to that due to molecular viscosity. However, in the regions 
away from the wall, they are comparable. Finally, in Fig. 18^ the cascade 
term, ” < > , is also plotted. Near the wall, its magnltuJe is 

larger than the dissipation due to eddy viscosity, and away from the wall they 

are of the same order of magnitude. Therefore, as was pointed out in Sect. 2, 

— ” 3 

inclusion of < u t ^ 1° modeling assumption for the subgrid scale 

1 dXj ij 

stresses is not recommended. As long as this term can be evaluated explic- 
itly, one should do so. 


7 . Detailed Flow Structures 

In this section we shall Investigate the detailed structure of the com- 
puted flow field. This will be done by examining contour plots of instanta- 
neous velocity, pressure, and vorticity field, and by tracking passive 
particles in the flow. The latter approach is a simulation of laboratory 
flow-visualization experiments using hydrogen-bubble wire. 

Figure 19 shows the contour plot of u in the x-z plane at y^ * 
6.26 and at the non-dimensional time, t • 4.3. In all the contour plots 
shown here, positive values are contoured by solid lines and negative values 
are contoured by dashed lines. In addition, all the instantaneous plots are 
obtained at the non-dimensional time t * 4.3. The distinct feature of the 
flow patterns in Fig. 19 is the existence of highly elongated regions of high- 
speed fluid (u > 0) located adjacent to the low-speed regions. This pic- 
ture of the flow in the vicinity of the wall ia in agreement with laboratory 
observations* In their visual studies, Runstadler et al. (1963) and, more 
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recently, other investigators have clearly doionstrated that the viscous 
sublayer consists of coherent structures of high- and loir-speed streaks alter-* 
nating in the spanwise direction. These studies have also shown that the 
streaks are the unique characteristic of the wall- layer turbulence, and they 
are absent in the regions away from the walls. Figure 20 shows the contour 
plot of u in an x-z plane far away from the wall (y'** " 400). In this 
region, in agreement with experimental observations, it is clear that the 
screaks and, for Chat matter, any definite organized structures are absent in 
the computed flow patterns. 

In Fig. 19, one can distinguish several localized regions of very high- 
speed fluid (large concentration of solid lines) that are located on the high- 
speed streaks. Figure 21 shows the corresponding contour plot of pressure 
fluctuations, obtained at the same vertical location (y'*' > 6.26). It can be 
seen that, in contrast to u , the pressure patterns are not elongated in the 
streamwlse direction. However, the regions of high-pressure fluctuations are 
generally located in the vicinity of the "pockets" (see Falco, 1978) of high- 
speed fluid. This correspondence together with examination of the contour 

plots of V (see below) suggest that these pockets are "quasi-stagnation" 

regions which are formed as a result of the arrival of high-speed fluid to the 
wall layer. Moreover, the contour plots of normal and spanwlse velocity 
fluctuations stiow that, like the pressure patterns, they do not exhibit 
elongated streaky structures. These observations imply that the wall layer 
may be viewed as a bed of low-speed fluid that is constantly subjected to the 
arrival of energetic eddies from the layers above. These energetic eddies 

(with the help of the strong mean shear) form =the high-speed streaks In the 
wall region. 
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Figure 22 shows the coutour plot of spanwise vortlclty fluctuation, 

n 

in the same x-s plane as in Fig. 19. Virtually all the regions with large 
vortlclty fluctuations are associated with negative (large concentration 
of dashed lines) . In these regions , the streamwise velocity profile has ex- 
cess momentim with respect to the mean velocity (i.e., u >0). It should 
be pointed out that, in the vicinity of the wall, the relatively large posi- 
tive values of S(u ) (Fig. 8) indicate that the existence of regions with 
u < 0 are more probable than those with u >0. however, the structures 
with large magnitudes of u are most likely associated with positive values 
of u . This is in agreement with the above observations. 

Figures 23 and 24 show the u and v patterns in an x-y plane (z “ 
4hj) which pass through the high-speed region in the lower left-hand corner 
of Fig. 19. In Fig. 24, a positive v (the solid lines) represents fluid 
moving in Che positive y-direction, and a negative v (the dashed lines) 
represents fluid moving in the negative y-dirsctiun. It can be seen that, in 
the vicinity of the walls, the high-speed fluid elements (u » 0) corres- 
pond to the sweep event, i.e., v < 0 near the lower wall, and v > U near 
the upper wall. On the other hand, the lo«r-speed fluid elements are generally 
being ejected from the wall regions. Clearly, both the sweep and ejection 
events have a positive contribution to the production of turbulent kinetic 
energy. One of the distinct features of Fig. 23 is that the high-speed struc- 
tures near the walls are Inclined at oblique angles with respect to the walls. 
This is the consequence of the action of mean shear on any fluid element from 
the outer layers that is moving toward the walls. Similar large-scale struc- 
tures have been identified in the laboratory by Rajagopalan and Antonia 


31 


(1979). In their neasurenants , they report the mean angle of inclination of 
theae structures to be 13* . 

Figures 25 and 26 ahow the contour plots of u and ^ in a y-a 
plane (x ■ 0). In these plots, only the lower half of the channel is shown. 
Throughout a significant portion of the region displayed, there is a negative 
correlation between u and v. Some of the regions where the correlation 
between u and v is negative extend from the wall region to the channel 
centerline. In the wall region, the vertical and spanwise extent of the 
eddies is significantly smaller than in the regions away from the wall. In 
particular, near the wall in Fig. 24 the array of hlghr and low-speed fluid is 
clearly discernible. 

Figure 27 shows the u , v , and w patterns in the close vicinity of 
the wall (y^ < 46, < 0.072) in the same y-z plane as in Figs. 25 and 
26. Here, the region near the wall is magnified, and hence the contour lines 
are highly distorted. In Fig. 27a, the mean spacing between two adjacent 
high-speed streaks (or low-speed ones) is about 250 in the wall units. The 
mean streak spacing can also be obtained from the ^^ll^^w “ 0.025, rj) pro- 
file in Fig. 7. In this figure, the negative peak occurs at r^ “ 125. This 
is the distance between two adjacent high* and lotf-speed streaks. Therefore, 
the corresponding distance between two high- (or low-) speed streaks is about 
250 in the wall units. These two values are, surprisingly, in good agreement 
with each other but are considerably larger than the generally accepted value 
of » 100. Therefore, as was pointed out in Sect. 3 for the Reynolds num- 
ber considered in this study, the computational grid resolution is inadequate 

to resolve the screaks at their proper scale. However, as we have seen, the 
computed flow patterns, in the wall region, do exhibit the streaky structures 
at the finest scale permitted by the grid. The k^(y^ - 0.025 ,t 3 ) profile 
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fron Case 1 in Table I shove that the awan streak spacing in that calculation 
is about 330 in the wan units* Thus» there is a definite iaprovaaent in the 
computed streak spacings with refinement of the computational grid resolution. 

In Fig* 27b, one can see an array of positive and negative regions of 
v-contour lines that correspond to fluid moving away from and toward the wall. 
Intense shear layers are located at the interface between the energetic fluid 
streams moving toward and away from the wall* These shear layers may undergo 
Helmholtz** type instabilities in the y-z plane that result in the formation 
of etreamwlse vortices. These vortices can clearly be identified in Fig* 28, 
where the contour plot of in the same y-z plane is shown. 

Comparison of Figs. 27b and 27c demonstrates that, in the close vicinity 
of the wall (y^ < 10), the high-speed vertical streams with negative normal 
component of velocity produce a flow pattern similar to that of a Jet impinge- 
ment on a plate. On the other hand, the high-speed vertical streams with 
positive normal component of velocity are formed from two streams with op- 
posite velocities in the spanwise direction. Since the high-speed fluid 
elements arriving at the wall region are more energetic than the viscous- 
dominated fluid moving away from the wall, there is a net transfer of energy 
from the normal component of turbulence intensity to the horizontal components 
(the splatting effect). This appears to be the reason for the behavior of the 
pressure-strain correlations in the vicinity of the wall (Fig. 11). In addi- 
tion, it should be noted that the impingement of fluid from outer layers on 
the wall leads to stretching of spanwise vorticity fluctuations (as well as 
streamwise vorticity) which can be an important mechanism for its amplifica- 
tion. 

Near the wall, in Fig. 27c, one can see large gradients of the spanwise 
velocity component in the normal direction* On the other hand, due to the 
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wall-suppression effect, the vertical component of turbulence intensity is 
degraded there. Large values of du/3y lead to high values of streamwise 
vorticity fluctuations in the neighborhood of the wall. In Fig. 28, it can 
be seen that the regions with large are concentrated near the wall. 

Here, two distinct areas can be identified: the first is slightly above the 
wall (10 < y*^ < AO) , where the large amplitudes of are due to revolving 

fluid elements Induced by the Intense shear layers shown in Fig. 27b^ the 
second Is in the iomediate neighborhood of the wall (y*** < 10) , where the 
splatting effect and no-slip boundary conditions lead to large values of 
3w/dy and consequently ( 0 ^^. In Fig. 9, the profile of < attained 

its maximum at the wall and displayed a local maximum at y*** » 30. The maxi- 
mum at the wall is a result of the splatting effect, and the local maximum is 
located in the first region described above. 

So far, we have examined the eddy structure of the turbulent channel 

flow by considering two-dimensional contour plots of instantaneous velocity, 
pressure, and vorticity field. To gain a better insight into the unsteady 

dynamics of the flow, a computer motion picture simulating flow-visualization 
experiments with hydrogen-bubble wires was made. Several sequences of film 
were generated. At regular intervals in each sequence “ 0.015 in non- 

dimensional time units) , 128 particles were generated along a line either 

parallel or normal to the walls. These particles were followed until the 
memory capacity of the graphic display unit was depleted. Here, we briefly 
discuss some of the still photographs taken from the film. 

Figure 29 shows tlie particles generated along a line parallel to the 
z-axis (“z-wire") and located near the lower wall (y*^ ■ 12). In this figure, 
the wall- layer streaks are clearly evident. Un several occasions when viewing 
the motion picture, it was observed that the particles generated near the wall 
were violently ejected to regions as far away from the wall as y*** ^ AOO. 
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Figure 30 ehowe the time history of the particles generated along a 
"s-iflre'* which Is located distant froa the wall (y'^ ■ 319). It can be seen 
that the coherent streaky structures that are the characteristic of wall layer 
turbulence are absent In the regions away froa the walls. 

In Fig. 31, the tine history of Che particles generated along a line 
normal to the walls Is shown. The fomatlon of Inflexional velocity profiles 
and strong shear layers near the walls Is very similar to the corresponding 
photographs obtained by Kim et al. (1971) in their flow-visualization 
studies. This resemblance Is even more pronounced In Fig. 32, where 128 
particles were generated along the same vertical line, as In Fig. 31, but 
extended from the lower wall to y ■ -0.5. Here, one can see several profiles 
with multiple inflexion points. In addition, In this figure, the formation of 
a streamwlse vortex with an axis of rotation which is tilted outward In the 
flow direction and its ultimate breakup are clearly discernible. 

8. Summary and Conclusions 

In this study, turbulent plane Polseullle flow has been numerically sim- 
ulated at a moderate Reynolds number. Host of the calculations were carried 
out with 316,096 grid points on the ILLIAC IV comuter. The agreement of the 
computed mean velocity profile and turbulence statistics with experimental 
data is good. 

The resolvable portion of the statistical correlations appearing In the 
Reynolds stress equations was calculated. The role and relative Importance of 
the various terms in these equations were discussed. 

The structure of the flow field was examined In sooe detail. It was 
found that. In agreement with experimental observations, the computed flow 
pattern In the wall region was characterized by coherent structures of low- 
and high-speed streaks alternating In the apanwlse direction. In this region. 
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Che large-anpiltude, screamwlse velocity fluctuacions were priaarlly due to 
Che arrival of hlglrapeed fluid eleaenca from adjacent layers. The regions 
with large-amplitude, screamwise vorclcity, were concentrated near the 

wall. Slightly above the wall, these regions contained revolving fluid ele- 
ments induced by strong shear layers in the cross-atream plane. In the imme- 
diate neighborhood of the wall, the splatting effect led to large magnitudes 
of and instigated transfer of energy from the normal component of 

turbulent kinetic energy to the horizontal components. 

With three-dimensional, time-dependent, numerical simulation of turbu- 
lence, one is capable of obtaining detailed, instantaneous information about 
the flow at thousands of spatial locations. This information can effectively 
be used to study the structure and statistical properties of the flow and 
their relation to each other. Furthermore, with the aid of computer graphics 
and the ability to move back in time and recreate an event in the flow aftei.* 
it has already been observed, one has the unique opportunity to study the 
mechanics of turbulent shear flows. Thus, with the anticipated advances in 
computer technology, it is expected that in the near future numerical simu- 
lation of turbulent flows will make Important contributions to turbulence 
research. 


This work was carried out in cooperation with the Thermo- and Gas-Oynamic 
Division of the Ames Research Center, NASA. We are Indebted to our 
colleagues, A. Leonard, R. S. Rogallo, and A. A. Wray, Anes Research Center, 
and J. H. Ferziger and W. C. Reynolds, Stanford University, for numerous 
helpful discuraions during the course of this study. 
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experimental data. 
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Fig. 3* Resolvable and total turbulence shear stress 
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Fi^« 4. Resolvable turbulence intensities and coaparison with experimei;tai 
data . 
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Fig. 5. Resolvable turbulence IntenelCiee in the vicinity of the lower wall 
and comparison with experiaentel date. 
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FlS> 13* Balance of the raaolvabla porcloa of the atraanwlaa coaiponaw 
turbulent kinetic energy ^ , production; Q , convection; 
caacade; V • velocity-preaaure gradient; q , diaaipation; < 
viacoua dlffuaion* 
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the resolvable turbulence kinetic 
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Fif. 21. Contours of P in the ic-z pleoe at 
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Pig. 22. Contours of spanvlse vortlclty fluctlon, u J , in the x-z plane 
At - 6.26. 
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Fig. 24. Contours of v In the x-y plane at r * 4hj . 
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Contours of u In the y-z plane at x ■ 0. (Only the lower half 
of the channel Is shown.) 


Fig. 26. Contours of v In the y-z plane at x ■ 0. (Only the lower lulf 
of the channel la shown.) 





Contours of u (a) 
♦6, y^6 < 0.072) 
















Fig. 29. Panicles generated from a "z-wirc" located at y* ■ 12. 
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ria» 31* Particles generated fron a vertical wire extended between the two 
channel walls; (a) two-dlnenalonal view; (b) thcee-dlwenslonal 
view. 
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